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SUMMARY 


The problem of extracting frequency, damping, amplitude, and phase infor- 
mation from unforced transient response data is considered. These quantities 
are obtained from the parameters determined by fitting the digitized time- 
history data in a least-squares sense with complex exponential functions. The 
highlights of the method are described and the results of several test cases 
are presented. The effects of noise are considered both by using analytical 
examples with random noise and by estimating the standard deviation of the 
parameters from maximum-likelihood theory. 


INTRODUCTION 


One of the fundamental tasks in flutter testing is the determination of the 
frequency and damping of aeroelastic modes. Transient response or free decay 
records are often used for extracting this information and may be generated 
directly by a method such as the resonant dwell and cut (e.g., see ref. 1), or 
indirectly through the use of autocorrelation or randomdec types of data- 
reduction techniques (refs. 2 and 3). Graphical or manual techniques have often 
been used to determine frequency and damping, but, with the widespread use of 
automated data-reduction procedures, numerical curve-fitting techniques of com- 
plex exponential functions or damped sine waves are frequently used. There may 
be strong interactions between the curve-fitting method and the data-collection 
process, especially in the areas of record length requirements and specifica- 
tions of noise level and distortion. Several procedures are currently available 
for the curve-fitting process (refs. 4 to 6). The purpose of this paper is to 
describe a method that takes a somewhat different approach from the previous 
works. The emphasis here is on developing a nearly real-time digital technique 
that is not only computationally fast but is also stable in the presence of 
real-world noise or contamination effects. A simple direct search technique 
for obtaining a least-squares fit using exponential functions has been developed 
and is presented. The application to several test cases is presented and dis- 
cussed. Some effects of measurement noise are evaluated by comparing test-case 
results for different signal-to-noise ratios, and by developing estimates of the 
standard deviations of the parameters from maximum-likelihood theory (ref. 7, 
e.g.) . 


It should be kept in mind that although in a practical engineering sense 
the use of exponential functions for the analysis of data may be satisfactory, 
the aeroelastic equations are not strictly constant-coefficient ordinary 
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differential equations (ref. 8) and may involve other functions. Furthermorej 
the extrapolation of damping measured at subcritical conditions to flutter may 
also have shortcomings. For example^ a case presented in reference 9 indicated 
a slope and curvature away from a flutter crossing in a plot of damping against 
velocity, even up to within 2 percent of the flutter speed. 


SYMBOLS 


coefficient in curve fit, the offset or static value (eq. (1)) 
coefficient of kth cosine term in curve fit (eq. (1)) 
coefficient of kth sine term in curve fit (eq. (1)) 
mean-squared error (eq . (2) ) 
expected value (eq. (3)) 
frequency, Hz 

ith data point of digitized time history 

data point index, 1 to N 

parameter index 

number of modes in curve fit 

modal index, 1 to K 

number of data points in digitized time history 

output error covariance matrix 

parameter sensitivity matrix 

time, seconds 

velocity 

flutter velocity 

curve-fitting expression (eq. (1)) 
fraction of critical damping 
damping coefficient (eq. (1)) 
f r equency , r ad / s ec 



ANALYSIS 


Least“Squares Fitting Procedure 

Given a free decay record containing the response of one or more vibration 
modes in the form of a digitized time history, the problem is to determine the 
modal damping, frequency, amplitude, and phase of each mode, A least-squares 
curve fit is made with complex exponential funetions (or damped sine waves) in 
the form 

K -n t 

Y(t) = aQ + I e ia^ cos sin Wj^t) (1) 

by minimizing the squared- error difference between the output fit Y(ti) and 
the input time history f^. The error is given by 

N 

E = I [Y(t ) - f r (2) 

i=l 


Inspection of equation (1) shows that if r|j^ and O)]^ are preassigned, it 
is possible to compute ag, and by solving a linear least-squares prob- 

lem. The nonlinear parameters r\y^ and 00 ]^^ must be determined by some type of 
search algorithm. Although this is a standard nonlinear, unconstrained optimi- 
zation problem for which several methods are available for trial, for simplicity 
a direct search technique is used to search the coordinate space until 

the values that minimize equation (2) are obtained. At each step, values for 
rik and oaic are determined, the small linear system solved, and the error 
recomputed. 

The technique has been programed for the Xerox Sigma 5 digital computer. 

In the program, the coordinate stepping process proceeds as follows: 

(1) A starting set of coordinates (k = 1, . . ., K) and a starting 

step size are furnished to the program. 

(2) The error E is computed at and at 4K additional points 

obtained by adding and subtracting the step size to or from each value of ri|^ 
and o)|^« If the central error E is less than any of the 4K peripheral 
values of E, the step size is reduced by 75 percent, and the calculations are 
repeated . 

(3) Otherwise, the point that gave the lowest value of E is taken to be 
the new central point, and the step size is increased by 10 percent. 

(4) The procedure is terminated when either the step size has been reduced 
below a preassigned threshold or a preassigned number of steps have been 
executed ^ 
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The method requires starting values for and oajr. For a single mode the 

starters can be arbitrary* However ^ for the multiple-mode cases the computer 
time can be significantly reduced by choosing good starters* The following 
procedure has been found to be a reasonable way of getting starters for 
multiple-mode cases: 

(a) Generate a one-mode solution using arbitrary starters. 

(b) Compute the difference between the one-mpde solution and the input 
data, that is, the output error. Then generate a one-mode fit to the error. 

(c) Use the and 0)|^ values from steps (a) and (b) as the starters 

for the two-mode solution. 

(d) For higher modes, steps (b) and (c) are repeated using the difference 
between the current multiple-mode solution and the original data to estimate 
the next higher mode. 

Although this procedure is cumbersome, it appears to be stable and, at 
least for the single-mode case, surprisingly fast. It would also be very help- 
ful to set the method up on an interactive basis similar to the technique 
described in reference 10. 

One of the schemes in the literature is referred to as Prony’s method 
(ref. 4). It computes and OO]^ by solving a 2K-order polynomial equation 

whose coefficients are determined from a least-squares process. The solution 
for the coefficients Sq, a^^, and is then determined by a linear least- 

squares procedure, as is done here. Since this method is elegant and computa- 
tionally efficient, it was examined during the present study. However, it has 
been the authors' experience that although Prony's method works well for perfect 
data, it is so sensitive to real-world noise that it is essentially useless even 
for generating starters for the search algorithm. 


Uncertainty Levels of Estimated Parameters 


The standard deviations of the estimated parameters, or uncertainty levels, 
can be determined from maximum-likelihood theory (ref. 7, e.g.). This type of 
estimate has provided some useful results in the field of stability and control 
(e.g., ref. 11). Assuming only measurement noise that is Gaussian and white, 
the expected variance of the parameter vector is 


E[A|J 







(3) 


where S is the parameter sensitivity matrix, is the output error covari- 
ance matrix, here a constant, and T denotes matrix transpose. The parameter 
vector p is made up of Bq, a^,, b|j^, and and the sensitivity matrix 

8Y(ti) 

is given by Sj (tj[^) = ^ — . These elements of the sensitivity matrix can be 
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calculated by directly differentiating equation (1) j and the variance is a 
normal output parameter, ThuSj for a single channel of data^ as considered 
here, these parameter uncertainty levels can be readily calculated after the 
curve-fitting process is completed. 


RESULTS AND DISCUSSION 


The curve-fitting method has been applied to three sets of data as test 
cases. The first case is a calculated damped sine wave with noise added with a 
random- number generator. The true answer is thus known. The second case is 
wind-tunnel data from the dynamic calibration of an aircraft gust vane. The 
third case is a set of data consisting of the subcritical randomdec signatures 
of the response to input noise of a two-dimensional flutter model that was 
implemented on an analog computer. 


Analytical Test Case 


The calculated data for the analytical test case with no added noise are 
shown in figure 1(a) and are compared with the fitted curve, which is exact in 
this case. For this case, the analytical input function was a single mode with 
offset and is given by 


Y(t) - 1 - cos 30t 

The curve fits for various levels of random noise are shown in figures 1(b) 
to 1(d). The noise level is def ined as the rms level of the Gaussian noise and 
is given as the fraction of the maximum amplitude of the mode that is 1. The 
results of the curve fit are summarized in figure 2. Only modest degradation 
of the results is shown for reasonable values of noise level of up to 0.10 or 
0.20. Also shown, as brackets on the points, are the standard deviations of 
the parameters, or uncertainty levels, calculated from equation (3) using the 
results output from the curve-fit procedure. In this case the exact modal 
parameters are known and it is possible to calculate a predicted uncertainty 
level from the exact parameters by assuming that the output error covariance is 
the value for the noise only. These predicted levels are shown as dashed lines. 
Both results give a good indication of the actual scatter, and thus the confi- 
dence level, with noise level. It might be noted that the effect of noise is 
larger on the coefficients aj and than on the damping, frequency, or 

offset. Thus, one must be more cautious in using the magnitude and phase 
information from such procedures. 

These results amply demonstrate that the algorithm works well in the pres- 
ence of random measurement noise. It has been the authors’ experience, however, 
that a test case of this type does not Indicate that a method will be satisfac- 
tory in practice. The noise here is random with zero mean, whereas in the real 
world, the effects of frequency drift, meandering means, and harmonic distortion 
are more severe. 


47 



Gust Vane Data 


The two cases considered are from a wind-tunnel dynamic calibration of a 
light balsa vane used to sense atmospheric turbulence on an aircraft. The vane 
was mechanically displaced and released repetitively. Since the response to 
tunnel turbulence was a sizable fraction of the total response and the release 
conditions somewhat ill defined, the transients were ensemble averaged. The 
background noise was thus diminished, but a pure step response was not obtained. 
The two cases considered are called the low-damping case and the high-damping 
case (although the low-damping case is relatively highly damped by structural 
standards) . The data for the low-damping case and the one-mode fit are shown 
in figure 3(a). The fit is reasonable, but there is some systematic deviation, 
particularly near the rightmost portion of the data. A two-mode fit was com- 
puted and, as shown in figure 3(b), gives a significantly Improved result. The 
results of the one- and two-mode fits are summarized as follows: 

One-mode fit 

Y(t) = 0.0444 

<3t 

-e (0.239 cos 140t - 0.452 sin 140t) 


Two-mode fit 

Y(t) = 0.0450 

-e-35'lt(o.i37 cos 133t - 0.372 sin 133t) 

-e-34-6t(o.086 cos 190t - 0.026 sin 190t) 

As compared with the one-mode results, the two-mode data indicate that the off- 
set is nearly the same, the frequency of the first mode is reduced by about 
5 percent, and the damping is reduced by about 10 percent, along with sizable 
changes in the coefficients of the first mode. The physical significance of the 
second mode is not clear in this case; it may be low-frequency noise that has 
not completely averaged out in the ensembling process. However, it is thought 
that the results for the first or principal mode obtained in the two-mode fit 
are more representative of the system response. 

The results for the highly damped case are presented in figures 4(a) and 
4(b). The results and trends are similar to those of the low-damping case. 

This case is a particularly difficult one to analyze, as it has high damping, a 
large offset, and a low-frequency distortion. The algorithm of this paper 
appears to give a reasonable result for this case. 


Randomdec-Analog Flutter Data 

Some subcritical randomdec signatures of the response of a two-dimensional, 
two-degree-of-freedom flutter model to input noise on an analog computer are 
also treated. The mathematical model and test setup were the same as those of 
the investigation of reference 9, The signatures and a one-mode curve fit are 
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shown in figures 5(a) to 5(d) as the velocity approaches flutter* As flutter 
is approached the signatures show little scatter or distortion^ in contrast to 
the lower velocities. The signatures contain two modes^ but the lower frequency 
mode is apparently unconverged in the randomdec procedure and could not be ade- 
quately resolved by the curve-fit procedure* The results for the flutter mode 
are compared with the exact solution in figure 6* The agreement is quite good 
in both frequency and dampings with the flutter speed underpredicted by less 
than 1 percent which is within the expected accuracy of the analog setup. 

Thus 3 the curve-fit procedure appears to be a practical means of analyzing 
randomdec signatures* 


CONCLUDING REMARKS 


A least-squares curve-fitting procedure to extract frequency, damping, 
amplitude, and phase information from free decay records has been presented. 

The method appears to be stable and to give reasonable results in the presence 
of noise. Some of the effects of noise on the parameter estimates can be 
assessed by calculating the uncertainty levels from maximum-likelihood theory. 
The method is relatively fast for a one-mode fit, generally requiring 5 to 
15 seconds on a Xerox Sigma 5 computer (which would be about 1 second on a 
CDC 6600 computer) and thus is a candidate for a real-time method. The two- 
mode solution, however, requires 2 to 5 minutes, and a three-mode solution is 
very long to calculate. Further work is needed to accelerate the multiple^mode 
calculations. It would also be very helpful to set the method up on an inter- 
active basis. Currently, the only multiple-channel capability is to fit each 
channel of data separately, determine a weighted mean for frequency and damping, 
and then recalculate the coefficients for each channel. This procedure may be 
satisfactory for engineering purposes, but the development to a true multiple- 
channel method may be desirable. 
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Figure 1.- Curve fit for analytical case 
with random noise. 
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Figure 2.- Results of analytical case with random noise. 
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Figure 3.- Gust-vane curve fit for 
low-damping case. 
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(a) One-Taode fit. 



Figure 4.- Gust-vane curve fit for 
high-damping case. 







